Three-body equations of motion in successive post-Newtonian approximations 
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There are periodic solutions to the equal-mass three-body (and TV-body) problem in Newtonian 
gravity. The figure-eight solution is one of them. In this paper, we discuss its solution in the first 
and second post-Newtonian approximations to General Relativity. To do so we derive the canonical 
equations of motion in the ADM gauge from the three-body Hamiltonian. We then integrate those 
equations numerically, showing that quantities such as the energy, linear and angular momenta 
are conserved down to numerical error. We also study the scaling of the initial parameters with 
the physical size of the triple system. In this way we can assess when general relativistic results 
are important and we determine that this occur for distances of the order of lOOM, with M the 
total mass of the system. For distances much closer than those, presumably the system would 
completely collapse due to gravitational radiation. This sets up a natural cut-off to Newtonian 
Ai'-body simulations. The method can also be used to dynamically provide initial parameters for 
subsequent full nonlinear numerical simulations. 

PACS numbers: 95.10.Ce, 95.30.Sf, 45.50.Pk, 04.25. Nx 



I. INTRODUCTION 

The closest star to the solar system, Alpha Centauri, is 
a triple system, so is Polaris and HD 188753. Triple stars 
and black holes are common in globular clusters (T], 01 , 
and galactic disks. Triple black hole mergers can be 
formed in galaxy merger [3] and a triple quasar, repre- 
senting a triple supermassive black hole system has been 
recently discovered 

Full numerical simulations of black holes made possible 
only in the last couple of years have already produced nu- 
merous astrophysically interesting results, among them, 
the orbital hangup and respect of the cosmic censorship 
hypothesis for spinning black holes f^, 0, 0| , precession 
and spin-flips J} , and the discovery [8| of large recoil ve- 
locities in highly-spinning black hole mergers up to 4, 000 
km/s [9*1 . 

The 2005 breakthroughs in Numerical Relativity [l^, 
[Til . [T^ , not only provided a solution to the long stand- 
ing two-body problem in General Relativity, but it also 
proved applicable to the black hole - neutron star bina- 
ries [l^ and recently to the three (and N) - black holes 
systems [l^ . 

In general, the solution of three-body problem in New- 
tonian gravity can be chaotic. There are however, pe- 
riodic orbits in the problem of three equal masses on a 
plane. One of the most surprising solution is a figure- 
eight orbit. The three bodies chase each other forever 
around a fixed eight-shaped curve. This was found first 
by Moore [IB] and discussed with the proof of the exis- 
tence in Ref. [31 . Heggie also estimates the proba- 
bility for such systems to occur in a galaxy. 

Because of effects such as the perihelion shift, it was 
unclear if the figure-eight orbits would exist in a low post- 
Newtonian expansion, even if it consist of only conserva- 
tive terms. Imai, Chiba and Asada succeeded in obtain- 
ing the figure-eight solution in a first post-Newtonian or- 



der approximation by finding the general relativistic cor- 
rections to the Newtonian initial conditions. In Ref. [l^ 
they also estimated the periodic gravitational waves from 
this system. 

In Ref. [l^ was used the Euler-Lagrange equations of 
motion in an approximation to first post-Newtonian or- 
der. In our paper we instead assume the Hamiltonian 
formulation to derive the equations of motion. We start 
from the Hamiltonian given in Ref. 2Ql (with typos cor- 
rected in our Appendix). We derive the equations of 
motion in this formalism, which are different from those 
used in Ref. [l^ and have the virtue of explicitly satis- 
fying the constants of motion of the problem, and thus 
being more amenable to numerical integration. 

The paper is organized as follows. In Section [TTl we 
summarize the equations of motion to be solved numeri- 
cally in order to obtain the figure-eight orbits. The start- 
ing point is the three-body Hamiltonian in the first post- 
Newtonian approximation. In Section [Till we discuss the 
initial conditions for the figure-eight solutions. We study 
the scaling relation between the orbital radius and the 
linear momenta. From this analysis, we can estimate 
when general relativistic effects are important. In Sec- 
tion IIVI we extend our calculation to the second post- 
Newtonian order and in Section |Vl we summarize the re- 
sults of this paper and discuss some remaining problems. 
The 2PN three-body Hamiltonian is explicitly given in 
the Appendix. Throughout this paper, we use units in 
which c = G — 1 . 



II. EQUATIONS OF MOTION 

As we mentioned in the introduction, the Newtonian 
configuration that leads to orbital braid figures can also 
be obtained within the Lagrangian approach, in the first 
post-Newtonian approximation, by finding the appropri- 
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ate corrections to the initial data [l9| . 

Here we will consider the Hamiltonian formulation 
since it generates equations of motion that conserve the 
energy, linear and angular momenta. This is crucial to 
reach high accuracy in the numerical integrations, which 
is needed to keep good track of the orbital motion, that 
in the three body problem might be chaotic. 

The Hamiltonian {H — Hn + HipN + H2pn) for the 
three body problem in the second post-Newtonian ap- 
proximation is given next. Note that since gravitational 
radiation only enters at 2.5PA'^-order and higher, the cur- 
rent analysis applies to conservative systems. 

The Newtonian Hamiltonian is given by 
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where a, b and c run over 1, 2 and 3. We have used the 
notations; x^f, = x^ - Xf,, rah = |xafc|, Uab = Xab/rab, 
Pa — Pa ' Pa ^nd thc dot (•) means the inner product. 
The Hamiltonian for the second post-Newtonian order is 
given in the Appendix. 

We then obtain the canonical equations 



iPa)i 



dH 

d{qa)'. 
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dH 

d{Pa)i 
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where i denotes x, y or z. 

Explicitly, the equation of motion for thc first post- 
Newtonian order, are given for thc particle 1 by 
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(5) 



where to obtain the equation of motion for the particle 2 
(and 3), we change the subscripts as {1 ^ 2, 2 — > 3, 3 ^ 
1} (and {1 ^ 3, 2 ^ 1, 3 ^ 2}), respectively. 

We solved the above equations numerically for three 
body problems using a 10 digits precision implemented 
in Maple 10 with typical runs of a few seconds on a Lap- 
top. Since we use the canonical momentum in the cal- 
culation, the Hamiltonian H, the total linear momentum 
P = ^ Pa and angular momentum L = ^ x^ x p^ are 
conserved quantities. These represent useful checks of 
the accuracy of the numerical runs. 



III. THE FIRST POST-NEWTONIAN 
CORRECTIONS 

In the Newtonian figure-eight motion can be 

obtained from the following initial conditions p^ . i.e., 
the positions 1 and linear momenta p: 



1 = 



(2:3, 2/3) 
Pn 



(xi, Vl) = {-X2, -2/2) 

(97.00, -24.31) , 
(0, 0) , 

{pl.,pl)^{-2pl^-2py) 
(-0.09324, -0.08647) . 



{-2pl, -2pl) 



(6) 
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Here, we set mi = TO2 = = m = 1. For the above 
initial condition, the total linear momentum and angular 
momentum are zero. 

At the IPN order, we also impose the total linear mo- 
mentum P = and the total angular momentum L = 0. 
By these conditions, we find that each linear momentum 
is given by the relations 



P3 



2pi 



-2p2 



(7) 



Therefore, when we give the positions of the three parti- 
cles, and it is necessary then only to search numerically 
for p3. In order to obtain, pa, we make some iterative 
computations until the figure-eight is reproduced for a 
few orbits. 

In Figures [T] and [5] we show the relative error of the 
Hamiltonian conservation: 



AH{t) = 



H{t) - H{0) 
H{0) 



Figure [T] is estimated by using the orbit calculated in 
[igj and we observe that they lead to violations of the 
order of 3 x 10"'^. While, Figure [2] is derived by using 
the canonical equations derived in our paper and they 
display errors of the order of 10~^, growing linearly in 
time due to the propagation of numerical errors triggered 
by initial roundoff. 

Next we will discuss the scaling behavior of p3 when we 
change the initial separation as 1 — > Al, and hence the size 
of the orbit. Note that pa — > A~^/^p3 in the Newtonian 
limit as can be easily derived from the Hamiltonian in 
Eq. Il]) or the equations of motion Q. 

In Table HI we summarize our numerical findings for 
the IPN initial conditions for A from 1 to 100. We note 



TABLE I: The initial conditions and inclination angle. 



A 




iP3)y 


e 


1.00 


-0.09811067089 


-0.09490870640 


0.01535863098 


2.00 


-0.06754964265 


-0.06392246619 


0.007238984240 


5.00 


-0.04209168100 


-0.03934705365 


0.002786451510 


10.00 


-0.02961805051 


-0.02758150399 


0.001351084509 


20.00 


-0.02089989478 


-0.01941808121 


0.0006871250545 


50.00 


-0.01319661317 


-0.01225031026 


0.0002447024114 


100.00 


-0.009328662000 


-0.008654573162 


0.0001269692928 



that p3 with A = 1 is different from the value which are 
derived from the initial velocity of [l^. The value 6 in 
the table is the inclination angle of the principal axes. 
The principal axes of the IPN figure-eight motion are 
not along the x and y axes [l^ . 

In Figure [3l the figure-eight rescaled orbits with A = 
1, 10 and 100 are shown. Here, in order to display the 
general relativistic effects, we have used the coordinates: 
{xa{t)/X, ya{t)/X). We have chosen here the x-axis as the 
principal axis. We observe that the superposition of the 
A = 10 and A = 100 is suggestive that at those scales the 
general relativistic effects are very small while for A < 1 
they are dominant, but remainder gauge effects may also 
mask this effect because the orbits are not gauge invari- 
ant. A cleaner analysis can be made directly looking at 
the initial linear momenta scaling. 

By using the results of the runs in Table |T1 we pro- 
pose a fitting formula for jpsj inspired again in the IPN 
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Hamiltonian or the equations of motion 



|P3|fit(A) 



0.01617387234 0.002042558971 0.0004169461512 



A 



A3 



(8) 




FIG. 3: Figure-eight motions. We show A = 1 (solid line) 
A = 10 (dashed line) and A — 100 (dotted line). 




FIG. 5: The relative error of the fitting, 
and the first post-Newtonian calculations: 

|P3|fit(A) - |P3|n(A) 

Ip3|n(A) 

is 7% for A = 1, 0.6% for A = 10 and 0.07% for A = 100. 



FIG. 4: A-|p3| relation with points obtained numerically. 



In Figure m we show the fitting function and in Figure [5] 
we display the relative error |P3| — |p3|fit/|P3| , consis- 
tent with the form of an error generated in the numerical 
calculation. 

Independently in the Newtonian calculations, the A- 
P3I relation can be obtained from the initial condition 
in Eqs. ([6]) as 



Ip3|n(a) 



0.1271642973 

A172 ■ 



(9) 



Note that relative difference jpaj between the Newtonian 



IV. SECOND POST-NEWTONIAN 
CORRECTIONS 

It is interesting to verify if this kind of orbits also exists 
in the second post-Newtonian approximation to General 
Relativity, since they incorporate further effects of the 
curvature, but yet not gravitational radiation. The cal- 
culations are done by using the same method as for the 
first post-Newtonian order. In Table [III we summarize 
the initial conditions for each A from 1 to 100. We show 
the numerical errors as measured through the Hamilto- 
nian non-conservation in Figure [6l 

We find that we can approximate jpaj by the fitting 
formula 
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|P3|fit(A) 



0.01617654493 0.002017242451 0.0002017242451 0.0001054698539 



A 



A2 



A3 



A4 



(10) 



TABLE 11: The initial conditions and inclination angle for the 
second post-Newtonian case. 



A 




(P3)y 


e 




1.00 


-0.09759146109 


-0.09386471063 


0.01335212441 


0.1 ■ 


2.00 


-0.06746813797 


-0.06375625776 


0.006775950067 




5.00 


-0.04208326266 


-0.03933131483 


0.002713363325 


o.os 


10.00 


-0.02961805051 


-0.02757874584 


0.001340868765 




20.00 


-0.02089780479 


-0.01941808121 


0.0006733410290 


50.00 


-0.01319661317 


-0.01225031026 


0.0002447024114 


o.os 


100.00 


-0.009328662000 


-0.008654573162 


0.0001269692928 






FIG. 6: The relative error of the Hamiltonian conservation 
for the second post-Newtonian order calculations. 



There is a significant difference between the coefficient of 
l/A^ in Eqs. ([8]) and (fTO| . This is due to second post- 
Newtonian corrections entering in this coefficient, as we 
can verify from the form of the Hamiltonian. 

In Figure [71 we show the fitting function while its rel- 
ative error is given in Figure [S] 

Finally, we summarize the results by showing the dif- 
ference between the Newtonian, first and second post- 
Newtonian results in Figure [51 The second post- 
Newtonian effect is small but clearly not negligible for 
A = 1. 



FIG. 7: A-|p3| relation for the second post-Newtonian case. 
The points are obtained numerically. 



V. DISCUSSION 

In this paper we have used the figure-eight orbits as 
a theoretical lab to test the properties of the low post- 
Newtonian expansions of General Relativity. We have 
found that those closed orbits exists for three (and pre- 
sumably N) bodies. We have provided an improved first- 
post-Newtonian order formalism for deriving the equa- 
tions of motion that satisfy the Hamiltonian (the lin- 
ear and angular momenta) constraint to round-off error. 
The subsequent numerical evolution is well behaved dur- 
ing for more than t ~ 10,000m. We have also extended 
this analysis to the 2PN corrections, still giving a con- 
servative system of equations. In the process of finding 
the figure-eight solutions by trial of different initial mo- 
menta we also showed (numerically) the stability of the 
orbit against small perturbations. 

This method is particularly useful to determine, dy- 
namically (as an alternative to determine them through 
families of initial data [2T| ) , initial orbital parameters for 
subsequent full numerical evolution [T^ . when the holes 
are close enough that general relativistic effects can no 
longer be ignored. Note that our method fully takes into 
account the three-body post-Newtonian interactions un- 
like other simulations that approximate the problem in 
successive two-body problems [22[. 
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FIG. 8: The relative fitting error for tlie second post- 
Newtonian case. 




FIG. 9: Comparison of figure-eight motions for A = 1. The 
solid, dotted and dashed lines show the 2PN, IPN and New- 
tonian results, respectively. 



It is interesting to note here that the scahng fits (fTO|) 
give a practical way to determine when relativistic or 
Newtonian approaches are appropriate. For A = 1 we 
have that the ratio of the first coefficient, 0.01617654493 
(Newtonian) to the second coefficient 0.002017242451 
first-post-Newtonian is nearly 0.12/A and the second co- 
efficient to the third one 0.0002463605227 (dominated 
by second-post-Newtonian) is also approximately 0.12/A. 
This indicates that post-Newtonian corrections are im- 
portant. For A = 1 the distance between the initial bod- 
ies is 200m, what indicates that for nearly 67M with 
M « 3m the total mass of the system has strong post- 
Newtonian effects. For A ^ 1 Newtonian gravity should 
describe the system accurately, while for A < 1 general 
relativistic effects should be very important, eventually 
leading to the total collapse of the system. It is inter- 
esting to remark here that most of the A^-body codes 
use some sort of regularization of the Newtonian gravity 
for very close encounters [1^, instead the natural way 
to regularize these close encounters [l3l is given by the 
General Theory of Relativity, and as we show here, the 
post-Newtonian corrections are already non-negligible at 
separations of the order of lOOM. In any case, for most 
of the astrophysical encounters this is way too short dis- 
tance, but it can obviously be reached in systems involv- 
ing black holes and neutron stars. 
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APPENDIX A: THE SECOND 
POST-NEWTONIAN THREE-BODY 
HAMILTONIAN 

In this appendix, we give explicitly the Hamiltonian for 
the three body problem at second post-Newtonian order 
in the ADM gauge since there are some typos in the sum- 
mation of 20] . The equations of motion used in our paper 
can be derived straightforwardly from this Hamiltonian, 
but are too cumbersome to write down here. 
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H2PN = T7T > , ^"a +777 
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pI 



-10 



pI (nafc • Pb) 



1 



marrib 



Tab 
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(Pa • Pb) {nab ■ Pa) {l^ab ' Pb) 
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9 2 
PaPb 



[Pa ■ Pa) 



, (riab-Pa) (nab-pb) 



mf:mt 
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^1 ^ m„m,m. r + 144-2 ^""^ ' P^^' - 50^^ + 17^^ 

1 ^ ("ab • Pa) (ngfe " Pb) , ^ ^ ("ab " Pb) ("ab ' Pc) , (ngb • Pfc) (llae ■ Pc) 

+ 14 + riab • Hoc f 

marrib mbrric mbUic ) 

mgmbmc ( {ngb ■ Po) (iiac ■ pc) ^ ^ {ugb ■ Pb) (ngc ■ Pc) 
8 rlf, I manic manic 

, ^ Pc ("ac • Pc)^ , . (nab ■ Pc) (Hoc " Pc) I 

J 

+ 9 7 — I — I — A2Kb + "ac)Kb + '^cb)i^zrTi- - loz-z- 
^ a,bMa,bi''<^b + rbc + rca) ^ ^'^^c m„mc 

^ r^ PaiPbj ^ PciPcj ^ PaiPaj \ 

mgmb ml mi, J 



a 

mgmbmc /g Pa ■ Pc - (ngb ■ Pg) (ngb ■ Pe) 
^ {rgb + rbc + rca)rgb\ mgmc 



4 i: 

a,b^a,Cy^a.b 



,U-/-U, ^l^-/-U, . {J 

g Pa • Pb - (ngb ■ Pa) (ngb ■ pfc) _ ^ pI - (ngfe ■ Pe)^ _ pI - (rigfc ■ pg)^ I 
1 mgmbmc 1 rngmbm"^ ^ 1 m^mb 

^g,b#g,c#b '^g,b#g,C7ia ''"^ ''"^ ^ g,b#g 

3 mgmbmc 3 mgmbmc ^ 3 m^m^ 

a,b^a,c^a a,b^a,c^a,b a^b^a 

1 ^ T m^JTlbnic f I ^ O O 2 2 r» 2 / \ 2 r-*^ 

~ fil ~3-^ 1 18r„;,rge - eOr^ftr;,^ - 24r„(,rgc(rg5 + rbc) + BOrgbrgcr;,^ + 56rlbrbc 

°* a,b#a,c#g,b ^a*- ^'"^ 

-72rgbrL + 35rt + Sr^J - i !I« . (Al) 

g.bj^g «b 
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